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Abstract 

The melting curve of aluminium has been determined from to ^150 GPa using first principles calcu- 
lations of the free energies of both the solid and liquid. The calculations are based on density functional 
theory within the generalised gradient approximation using ultrasoft Vanderbilt pseudopotentials. The 
free energy of the harmonic solid has been calculated within the quasiharmonic approximation using 
the small-displacement method; the free energy of the liquid and the anharmonic correction to the free 
energy of the solid have been calculated via thermodynamic integration from suitable reference systems, 
with thermal averages calculated using ab-initio molecular dynamics. The resulting melting curve is in 
good agreement with both static compression measurements and shock data. 

1 Introduction 

The determination of the melting curves of materials to very high pressures is of fundamental importance to 
our understanding of the properties of planetary interiors; however, obtaining such melting curves remains 
a major challenge to experimentalists and theorists alike. In particular, the melting behaviour of iron 
is of great interest to the Earth science community, since knowledge of this melt transition would help 
constrain the temperature at the inner core boundary (about 1200 km from the centre of the Earth) which 
is currently uncertain to within a few thousand degrees. Although several attempts have been made to 
obtain the melting curve of iron, experimentally and theoretically determined melting curves vary widely 
with significant disagreement between static compression measurements m, I, i, shock data |, § and first 
principles calculations |8|, ^, 10 1. Consequently, the true nature of the melting curve of iron remains 



in some dispute. 

In order to test the reliability of the theoretical techniques used in our previous work on iron and to 
validate further the reported melting curve [^, ^] , we have calculated the melting curve of aluminium, for 
which there is a plethora of ambient experimental data (e.g. ||ll|l), and for which the experimental melting 
curve has recently been measured [|l2|, p^ . 

In the past, a number of theoretical approaches have been used to investigate the melting behaviour 
of aluminium. Moriarty et al. [jl^ used the generalised pseudopotential theory (GPT) to calculate the 
free energy of both the solid and liquid. They treated the solid harmonically within the quasiharmonic 
approximation and for the liquid they used fluid variational theory, where an upper bound for the free 
energy is calculated from a reference system constructed within GPT. They obtained a melting curve 



to 200 GPa in fair agreement with more recently determined experiment data 13, 14], predicting a 
zero pressure melting temperature of 1050K compared to the experimental value of 933K |11]. Mei and 
Davenport [|l^ used the embedded atom model (EAM) based on an analytical potential fitted to the 
structural properties of aluminium. They calculated the free energies of the solid and liquid and obtained 



a melting temperature at zero pressure of 800K. Morris et al. [17| employed the same EAM model but 



they used phase coexistence to determine the melting temperature as a function of pressure, with results 
considerably lower than previous theoretical and experimental estimates; they obtained a zero pressure 
melting temperature of ~720 K. Straub et al. |jl^ used first principles calculations to construct an optimal 
classical potential, and used this potential to calculate the free energies of the solid and the liquid using 
molecular dynamics; they obtained a zero pressure melting temperature of 955 K. 

The first fully ab-initio determination of aluminium melting behaviour is that of de Wijs et al. p^ , 
who obtained the zero pressure melting point by calculating the free energy of the solid and the liquid 
entirely from first principles. Their calculations were based on density functional theory (DFT) using 
the local-density approximation (LDA) for the exchange-correlation energy. The free energy of the solid 
was obtained as the sum of the free energy of the harmonic solid, within the quasiharmonic approximation. 



and the full anharmonic contribution, calculated using thermodynamic integration |21| using the harmonic 
solid as the reference system. For the liquid they used thermodynamic integration with a Lennard-Jones 
fluid as the reference system. They obtained a melting temperature of 890 K. More recently, Jesson 



and Madden [22| used the orbital-free (OF) variant of ab-initio molecular dynamics and thermodynamic 



integration to calculate the free energy of liquid and solid aluminium. They found a melting temperature 



of 615 K, attributing the discrepancy with the DFT-LDA value of de Wijs et al. [19| to either the OF 
approximation or the pseudopotential used. 

In this paper we present the first fully ab-initio calculations of the entire melting curve of aluminium 
from to 150 GPa. Our calculations are similar in the general principles to those of de Wijs et al. [^] 
in the sense that we calculate the ab-initio free energies of both liquid and solid using thermodynamic 



integration, although we use the generalised gradient approximation (GGA) ||2^, 24] for the exchange- 
correlation energy. In addition to extending the calculations to a wide range of pressures, we also present a 
more efficient approach to the thermodynamic integration scheme, in which additional intermediate steps 
are introduced in order to minimise the computational effort. Finally, we discuss some possible limitations 
of the GGA. 

The paper is organised as follows: in section ^ we describe the ab-initio simulation techniques and the 
strategy to calculate the melting curve; in section ^ and ^ we describe the calculations of the free energy 
of the liquid and the solid respectively, and in section H we present the melting properties of aluminium. 

2 Ab-initio simulation techniques and strategy for melting 

In the present work, the aluminium system was represented by a collection of Al^+ ions and 3iV electrons, 
where N is the number of atoms. The ions were treated as classical particles, and their motion was 
adiabatically decoupled from that of the electrons via the Born-Oppenheimer approximation. For each 
position of the ions, the electronic problem was solved within the framework of DFT [pO| using the GGA 



of Perdew and Wang 24 1. Thermal electronic excitations were included using the standard methods 



of finite-temperature DFT developed by Mermin |25, p3, 27 1. The present calculations were performed 



with the code VASP |28[] which is exceptionally efficient for metals. The interaction between electrons and 
nuclei was described with the ultrasoft pseudopotential (USPP) method [29|. We used plane- waves with a 



cut-off of 130 eV. The Brillouin-zone was sampled using Monkhorst-Pack special points ||30[ (the detailed 
form of sampling will be noted where appropriate). The extrapolation of the charge density from one step 
to the next in the ab-initio molecular dynamics (AIMD) simulations was performed using the technique 



described by Alfe |31], which improves the efficiency of the calculations by almost a factor of two. The 
time step used in our simulations was 1 femto-second. 

To calculate the melting temperature we calculated the Gibbs free energy of both the solid and the 



2 



liquid as a function of pressure and temperature, Gs(-P, T) and Gi{P,T) and, at each chosen P, obtained 
the melting temperature, T^, from Gs{P,Tjn) = Gi(P,Tm)- In fact, we calculated the Helmholtz free 
energy F(y, T) as a function of volume and temperature, and the Gibbs free energy was obtained from the 
usual expression G = F + PV, where P = —{dF/dV)T is the pressure. The main problem in determining 
melting curves with this technique is the high precision with which the free energies need to be calculated. 
This is because the Gibbs free energy of the liquid crosses the Gibbs free energy of the solid at a shallow 
angle, the difference in the slopes being the entropy change on melting. For aluminium this is about 1.4 
fee/atom at zero pressure, which means that an error of 0.01 eV/atom in either Gg or Gi results in an error 
of ~ 80 K in the melting temperature. Therefore, it is important to reduce non-cancelling errors between 
the liquid and the solid to an absolute minimum. In the next sections we give a detailed discussion of the 
techniques that we have used to calculate the free energies of the liquid and the solid, and report what the 
controllable errors are: those due to k-point sampling, finite size, and statistical sampling. We also try to 
give an estimate of what the uncontrollable errors due to DFT-GGA may be. 



3 Free energy of the liquid 

The Helmholtz free energy F of a classical system containing A'" particles is: 

P = -k^T\ny^ |^dRi...dR^e-/^^(^-^-^)}, (1) 

where A = h/ {2711^1 k^TY I is the thermal wavelength, with M the nuclear mass, h the Plank's constant, 
k-Q the Boltzmann constant and (i = l/k^T. The multidimensional integral extends over the total volume 
of the system V. 

A direct calculation of F using the equation above is impossible, since it would involve knowledge of 
the potential energy C/(Ri, . . . RAr;T) for all possible positions of the N atoms in the system. We have 
used instead the technique known as thermodynamic integration [p^], as developed in earlier papers p^. 



33| , |19| , |34| . This is a general scheme to compute the free energy difference F — Fq between two systems 
whose potential energies are U and Uq respectively. In what follows we will assume that F is the unknown 
free energy of the ab-initio system and Fq is the known free energy of a reference system. The free energy 
difference F — Fq is the reversible work done when the potential energy function Uq is continuously and 
reversibly switched to U. To do this switching, a continuously variable energy function Ux is defined such 
that for X = 0,Ux = Uq and for X = l,Ux = U. We also require Ux to be differentiable with respect to A 
for < A < 1. A convenient form is: 

C/a = (1-/(A))C/o + /(A)C/, (2) 

where /(A) is an arbitrary continuous and differentiable function of A with the property /(O) = and 
/(I) = 1. The Helmholtz free energy of this hybrid system is: 

i^A = -feBrin|^ |dRi...dR^e~^^^(^--^-^)}, (3) 

Differentiating this with respect to A gives: 

dFx^_. jpUlydK,... dn^e-^^^(^^-^--^^\-P^) ^ /dUx\ 
dX ^ ^ Jy dRi . . . (iR^e-/5^A(Ri,...R.v;T) \ 5A /, ' ^ ^ 
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so 

1 



AF = F-Fo= I dX{^j^. (5) 







For our calculations we defined Ux thus: 

Ux = {l- A) Uo + XU. (6) 

Differentiating Ua with respect to A and substituting into Equation |5| yields: 

1 

AF = y' dX{U-Uo)x- (7) 



Under the ergodicity hypothesis, thermal averages are equivalent to time averages, so we calculated {■)^ 
using AIMD, taking averages over time, with the evolution of the system determined by the potential 



energy function U\. The temperature was controlled using a Nose thermostat |35, 36 1. It is important to 
stress that the choice of the reference system does not affect the final answer for F, although it does affect 
the efficiency of the calculations. The latter can be understood by analysing the quantity {U — C/o)a- If 
this difference has large fluctuations then one would need very long simulations to calculate the average 
value to a sufficient statistical accuracy. Moreover, for an unwise choice of Uq the quantity {U — Uo)x may 
strongly depend on A so that one would need a large number of calculations at different A's in order to 
compute the integral in Eq. |^ with sufficient accuracy. It is crucial, therefore, to find a good reference 
system, where "good" means a system for which the fluctuations of U — Uq are as small as possible. In 
fact, if the fluctuations are small enough, we can simply write -F — -Fq — (f^ ~ t^o)oi with the average taken 
in the reference ensemble. If this is not good enough, the next approximation is readily shown to be 0: 

F-Fo^{U- Uo)o - ^ {[U-Uo-{U- Uo)of)^ ■ (8) 

This form is particularly convenient since one only needs to sample the phase space with the reference 
system, and perform a number of ab-initio calculations on statistically independent configurations extracted 
from a long classical simulation. 

To evaluate the integral in Equation |^ one can calculate the integrand {U — Uo)x at a sufficient number 
of A and calculate the integral numerically. 



Alternatively, one can adopt the dynamical method described by Watanabe and Reinhardt |37]. In this 
approach the parameter A depends on time, and is slowly (adiabatically) switched from to 1 during a 
single simulation. The switching rate has to be slow enough so that the system remains in thermodynamic 
equilibrium, and adiabatically transforms from the reference to the ab-initio system. The change in free 
energy is then given by: 

T • 

AF= I dt^{U-Uo), (9) 



where Tgim is the total simulation time, X{t) is an arbitrary function of t with the property of being 
continuous and differentiable for < t < 1, A(0) = and A(Tsim) = 1- 

When using this second method, it is important to ensure that the switching is adiabatic, i.e. that Tgim 
is sufficiently large. This can be achieved by changing A from to 1 in the first half of the simulation, and 
then from 1 back to in the second half of the simulation, evaluating AF in each case; the average of the 
two values is then taken as the best estimate for AF, and the difference is a measure of the non-adiabaticty. 
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If this difference is less than the desired statistical uncertainty, one can be confident that the simulation 
time is sufficiently long. 

In our calculations we chose a total simulation time of sufficient length such that the difference in AF 
between the two calculations was less than a few meV/atom. We return later to estimate the errors in our 



calculations in Section 3.5 



As pointed out by Jesson and Madden [g^], a possible problem in the calculation of the thermodynamic 
integral is that the system Ux may be in the solid region of the phase diagram, even though the two end 
members Uq and U are in the liquid region. If this happens, the system can freeze during the switching, 
and the integration path is not reversible, leading to an incorrect result. For small systems the situation 
is even more problematic, since the phase diagram is not defined by sharp boundaries, and the system can 
freeze even if it is above the melting temperature of the corresponding system in the thermodynamic limit. 
We have ourselves experienced freezing of the system for some simulations at temperatures very close to 
the melting point; in order to avoid including the results from these simulations, we carefully monitored 
the mean square displacement and the structure factor of the system, and included only those simulations 
in which these two quantities clearly indicated liquid behaviour throughout the whole simulation. 

3.1 The reference system 

We mentioned earlier that the efficiency of the calculations is entirely determined by the quality of the 
reference system, i.e. by the strength of the fluctuations of AU = U — Uq. The key to the success of 
these simulations, therefore, is being able to find a reference system such that the fluctuations in AU are 
as small as possible. Based on the experience of previous work on liquid Al ||l^ and liquid Fe |^ we 
experimented with the Lennard- Jones (LJ) system and an inverse power potential (IP). Analysis of the 
fluctuations in AU indicated that the system which best represented the liquid was the IP: 

UiP = lj2^i\^i-^j\)' (10) 

where 

^(^) = ^- (11) 
The potential parameters B and a were chosen by minimising the quantity: 

[Uip{B,a)-U-{Uip{B,a)-U)f) (12) 

with respect to B and q, where () means the thermal average in the ensemble generated by the ab-initio 
potential. To investigate whether the optimum values for the potential parameters depended strongly on 
thermodynamic state, we performed the optimisation at the three thermodynamic states of the extremes 
of high P/T and low P/T, and also a point in between; we found that the single choice oi B = 246.67 and 
a = 6.7 (units of eV and A) was equally good for all states and we therefore used these two parameters 
for all our calculations. 

It may be surprising that such a simple inverse power potential can reproduce the energetics of the 
liquid with sufficient accuracy, since simple repulsive potentials cannot describe metallic bonding. One 
may think that a more realistic potential such as those based on the EAM |16, 35, 39, ^] would be more 



appropriate, since these potentials explicitly contain a repulsive and a bonding term. However, in our 
recent work on iron [Q] we tested the use of an EAM potential as a reference system and found that the 
bonding term is almost independent of the positions of the atoms, depending only on the volume and 
temperature of the system, and the fluctuations of the energy are almost entirely due to the repulsive 
term. Since the only relevance in this work is the strength of the fluctuations (Eq. [l^ ), little is gained by 
using an EAM rather then a much simpler inverse power potential. 
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3.2 Free energy of the reference system 

Consider the excess free energy of the IP, Fjp = Fip — FpQ, where FpQ is the Helmholtz free energy of 
the perfect gas and Fip the total Helmholtz free energy of the IP system. The very simple functional form 
of Uip makes it easy to show that the adimensional quantity Ffp/kBT can only depend non-trivially on a 
single thermodynamic variable, rather then separately on V and T: 

Ff^/kBT = /(C) (13) 

with 

C = B/V^'^kBT. (14) 



The free energy of the IP has been studied extensively in the past |41], but only for special values of the 
exponent a, which did not include our own a = 6.7. We have therefore explicitly calculated the free energy 
of our inverse power potential using thermodynamic integration as before, but this time we started from 
a system of known free energy, the Lennard-Jones liquid, whose potential function is given by 



The free energy of the Lennard-Jones liquid, Flj, has been accurately tabulated by Johnson et al. p^ . 
To calculate Fip — Flj = AFlj^ip we used simulation cells containing 512 atoms with periodic boundary 
conditions and a simulation time Tsim = 200 ps. We performed the calculations for C ranging from 2.5 
to 6.25, with steps of 0.25. The calculations were done at a fixed volume of 14 A^/atom and varing 



temperatures according to Equation 14. We carefully checked that the results were converged to better 
than 1 meV/atom with respect to the size of the simulation cell and the length of the simulations. To 
avoid truncating the inverse power potential at a finite distance we used the Ewald technique. Our results 
were fitted to a third order polynomial in C: 

/(C) = Ec.c. (16) 



The coefficients are: cq = 2.4333; ci = 27.805; C2 = — 5.0704;c3 = 1.5177, and the fitting function repro- 
duced the calculated data such that the errors in Fip were generally less than 1 meV per atom. 

As an additional check on the calculated free energy, we repeated most of the simulations using the 
perfect gas as the reference system, thereby avoiding the inclusion of any possible errors that may exist in 



the free energy of the LJ system reported in the literature |42]. For these calculations we used a different 
form for U\, namely: 

Ux = X^Uip (17) 
(the potential energy of the perfect gas is zero, so does not appear in the formula). So Eq. ^ becomes 

FiP - FpG = t '^HUip)x. (18) 





The advantage of using this different functional form for Ux is that the value of the integrand does not need 
to be computed for A = 0, where the dynamics of the system is determined by the perfect gas potential. In 
this case, since there are no forces in the system there is nothing stopping the atoms from overlapping, and 
the potential energy Uip diverges. Not computing the integrand at A = partially solves this problem, 
but for small values of A where the forces on the atoms are small, the atoms can come close together 
and the potential energy, Uip, fluctuates violently. However, we found that by performing long enough 
simulations, typically 1 ns, we could calculate the integral with an accuracy of ~ 1 meV/atom, and, within 
the statistical accuracy, we found the same results as those obtained using the LJ reference system. 
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3.3 Free energy of the ab-initio system 



To calculate the full ab-initio free energy of the liquid, Fuq, we used thermodynamic integration, starting 
from the IP system. The calculations were performed at 18 different thermodynamic states over a range 
of volumes (9.5-19.5 A^/atom) and temperatures (800-6000 K). To address the issue of k-point sampling 
and cell size errors in the free energy difference F^^ — Fip , tests were carried out on cells containing up to 
512 atoms and a 5 x 5 x 5 k-point grid (calculations on the largest 512 atoms cell were only performed 
with T-point sampling), at V = 16.5 and T = lOOOX. The free energy difference F^^ — Fip was calculated 
using the perturbational approach (Eq. ^), with sets of configurations generated using the IP potential. 
We found that a 64-atom cell with a 3 x 3 x 3 k-point grid was sufficient to get convergence to within 2 
meV/atom. However, we were reluctant to perform simulations using the desired 3x3x3 k-point grid 
(14 points in the Brillouin Zone (BZ)) since these calculations are extremely expensive. We found it more 
efficient to add one further step to our thermodynamic integration scheme: 



AFr^333 = i^333 -Fr = J dX (C/333 - Ur),, , (19) 

1 

where C/333 ^-nd Up are the ab-initio total energies calculated using the 3x3x3 k-point grid and F-point 
sampling respectively, and F333 and Fp are the corresponding free energies. To evaluate the free energy 
difference AFr_+333 we noticed that the difference C/333 ~ Up did not depend significantly on the position 
of the atoms, so the integral in Eq. ^ could be evaluated using the second order formula (Eq. |8|). Using 
a long T-point ab-initio simulation, we extracted up to 25 statistically independent configurations and 
calculated the ab-initio energies using the 3x3x3 k-point grid. To test this, we performed spot checks at 
two thermodynamic states, where we calculated the full thermodynamic integral -F333 — Fip using adiabatic 
switching with a switching time of ~ 2 ps, and found the same results to within a few meV/atom. 

The free energy difference AFjp^r = Fp — Fip was obtained by full thermodynamic integration between 
the ab initio and reference system using adiabatic switching (Eq. ^) with a switching time of 5 ps, which 
resulted in errors of 1(4) meV/atom in the low(high) P/T region. To test this, we also calculated this free 
energy difference at several state points by numerical evaluation of the thermodynamic integral (Eq. ^), 
with A = 0, 0.5 and 1; we found that this gave the same numerical answer to within our statistical errors. 

In summary, the free energy of the liquid was obtained from a series of thermodynamic integration 
calculations: 

i^liq = i^333 = i^LJ + AFlj_ip + AFjp^r + AFr_^333- (20) 



3.4 Representation of the free energy of the hquid 

The results of the calculations described in the previous section were fitted to a suitable function of T and 
V. In order to do that efficiently we expressed the free energy in the following way: 

i^iiq = i^iP +AF = Fip + AU' + {AF - AU') (21) 

where AU^ = C/* — Ufp, with the zero temperature ab-initio (free) energy of the face-centred-crystal 
(fee) and Ufp the inverse power energy. can be calculated very accurately, details of which will be given 
below in Section 4.1; Ufp has no errors. The remaining quantity AF — AU^ is a weak function of V and 
T, and was fitted to a polynomial in V and T: 

AF-AU' = J2(Y. «y (22) 

j=0 \i=0 / 

The fitting reproduced the calculated data to within 2 meV/atom. 
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3.5 Error estimates for Fuq 

The errors on Fip and AV^ are each less than 1 meV/atom (see Section [4.1| below). The part of the free 
energy that carries the largest errors is AF — AC/*, which we estimate to be 2(5) meV/atom at low(high) 
P/T. 



4 Free energy of the solid 

The free energy of the solid can be represented as the sum of two contributions: the free energy of the 
perfect non-vibrating fee crystal and that arising from atomic vibrations above zero Kelvin: 

i^sol = Fperf + Fvib. (23) 

The contribution to the free energy due to the vibrations of the atoms may be written: 

arm ~l~ -^anharm (^4) 

where i*harm is the free energy of the high temperature crystal in the harmonic approximation and -Fanharm 
is the anharmonic contribution. 



4.1 Free energy of the perfect crystal 



The free energy of the perfect crystal, -Fperf, was calculated as a function of volume and temperature. 
Calculations were performed on a fee cell at a series of volumes (9.5-19.5 A'^/atom representing compression 
up to ~ 150 GPa) and temperatures (up to 6000K) with a 24x24x24 k-point grid (equivalent to 1300 points 
in the irreducible wedge of the Brilloiun zone (IBZ)), which ensures convergence of the (free) energies to 
better than 1 meV/atom. At each different temperature we calculated the ab-initio (free) energy as a 
function of volume, and then performed a least-square fit of the results to a third-order Birch-Murnaghan 
equation of state: 



E{V) = Eo + -VqK 
e = ^(4-i^')- 



4/3 



V 



2/3 



+ 2 ^+2 



(25) 



The parameters £^0) ^o^Kq, and K' were fitted to a fourth order polynomial as function of temperature: 

4 4 4 4 

Eo{T) = Y,eo,ir-, V^{T)=Y,v^,,r- Ko{T) = Y^k^^,T- K' {T) =Y^k'^^X . (26) 

1=0 i=0 j=0 i=0 

The fitting reproduced the calculated energies to better than 1 meV/atom in the whole P/T range. 

4.2 Free energy of the harmonic crystal 

The free energy of the harmonic crystal is given by: 



' BZ ^ 



1 

24 



+ . . . dq 



(27) 
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where ujq^i{V,T) are the phonon frequencies of branch i and wavevector q, JIbz is the volume of the 
Brillouin zone, Ni is the total number of phonon branches and the dependence on temperature of tUq^j is 
due to electronic excitations. We truncate the summation after the first term, which is the classical limit 
of the free energy: 

* BZ ^ -1' / 

This is a justifiable approximation to make for two reasons: (i) the error in making such a truncation is 
very small (<1 meV/atom), and (ii) neglecting the higher order terms, i.e., the quantum corrections, is 
consistent with the liquid calculations where the motions of the atoms were treated classically. 

It is useful to express the harmonic free energy in terms of the geometric average (D of the phonon 
frequencies, defined as: 

lnO = j^Y.H^ni), (29) 

where we have replaced the integral I with the summation J2- This allows us to write: 

BZ q 

i^harm = 3A;Brin(/3to). (30) 



To calculate the vibrational frequencies Wq^j, we used our own implementation of the small dis- 
placement method ^. 

The central quantity in the calculation of the phonon frequencies is the force-constant matrix ^isa,jti3j 
since the frequencies at wavevector q are the eigenvalues of the dynamical matrix 0^^,1/3^ defined as: 

Dsa,tf3{(l) = ^jj^jyj^ ^isa,jtl3 exp [iq • (R° + n - R° - T^)] . (31) 

where is a vector of the lattice connecting different primitive cells, is the position of the atom s 
in the primitive cell and Ms its mass. If we have the complete force-constant matrix, then Dsa,ti3, and 
hence the frequencies lo^i, can be obtained at any q. In principle, the elements of ^isajt/s are non-zero 
for arbitrarily large separations | R^ + Tt — R^ — | , but in practice they decay rapidly with separation, 
so a key issue in achieving our target precision is the cut-off distance beyond which the elements can be 
neglected. 

In the harmonic approximation the a Cartesian component of the force exerted on the atom at position 
R^ + Ts is given by: 

Fisa = — ^ ^isa,jtf3 Ujtp (32) 

where Ujsp is the displacement of the atom in R^ + rj along the direction /?. The force constant matrix 
can be calculated via: 

^isaJtfS — {^^) 

Ujt(3 

where all the atoms of the lattice are displaced one at a time along the three Cartesian components by 
Ujtfs, and the forces Fisajt/s induced on the atoms in R? + are calculated. Since the crystal is invariant 
under translations of any lattice vector, it is only necessary to displace the atoms in one primitive cell 
and calculate the forces induced on all the other atoms of the crystal, so that we can simply put j = 0. 
The fee crystal has only one atom in the primitive cell, so only three displacements are needed. However, 
a displacement along the x direction is equivalent by symmetry to a displacement along the y or the z 
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direction, and therefore only one displacement along an arbitrary direction is needed. It is convenient to 
displace the atom along a direction of high symmetry, so that the supercell has the maximum possible 
number of symmetry operations. These can be used to reduce the number of k-points in the IBZ, minimising 
the computational effort. For an fee crystal this is achieved by displacing the atom along the diagonal of 
the cube. 

Tests for cell-size (64-512 atoms), displacement length (0.01-0.0005 fraction of nearest neighbours dis- 
tance) and k-point grid (up to 9 x 9 x 9) were performed at the two extremes of high P/T and low P/T 
state points. Convergence of the free energy to within less than 3 meV/atom was achieved using a 64-atom 
cell with a 0.001 fractional displacement and a 9 x 9 x 9 k-point grid (equivalent to 85 points in the IBZ 
of the supercell). Calculations were performed for y = 9.5 - 18.5 A and T = 500 - 6000i^, and ln(a)) has 
been fitted to the following polynomial in V and T: 

H^)=jl{jl<^^^y']TK (34) 

j=0 \i=0 / 

The fitting reproduced the calculated data within « 1 meV/atom. 



4.3 Anharmonicity 

To obtain the anharmonic contribution to the free energy of the solid we have again used thermodynamic 
integration. In this case a natural choice for the reference system could be the harmonic solid [|19|, but 
unfortunately this does not reproduce the ah initio anharmonic system with sufficient accuracy. A much 
better reference system is a linear combination of the harmonic ab-initio and the same IP used for the 
liquid calculations 0: 

Urei = aUlp + 6[/harm, (35) 

where the harmonic potential energy is: 

t^arm = 2 5Z Uisa^isa,jt/3 'UjiPi (36) 

and where Ujsp is the displacement of the atom in -|- along the direction (3, and ^isa,jtp is the force 
constant matrix. The parameters a and h are determined by minimising the fluctuations in the energy 
differences U-ref — U on a set of statistically independent configurations generated with C/ref- However, when 
we start our optimisation procedure we do not know C/ref , so we cannot use it to generate the configurations. 
We could use the ab-initio potential, but this would involve very expensive calculations. We used instead 
an iterative procedure, like in our previous work on iron ||^]. We generated a set of configurations using the 
harmonic potential C/harm and calculated the ab-initio energies. By minimising the fiuctuations of U^ei — U 
we found a first estimate for a and b, and we constructed a first estimate of Uref- We generated a second 
set of configurations using this U^cf, calculated the ab-initio energies and minimised again the fluctuations 
of f/rcf — U with respect to a and b. This procedure could be continued until the values of a and b no 
longer changed, but in practice we stopped after the second step, and found a = 0.95 and b = 0.12. We did 
not use the extra freedom in the choice of the inverse power parameters since we found that this reference 
system already described the energetics of the solid very accurately. 

The calculation of the anharmonic part of the free energy required, once more, two thermodynamic 
integration steps. In the flrst step we calculated the free energy difference F^f — -Fharm- These are 
cheap calculations since they involve only the classical potentials Ujp and C/harm! the simulations were 
performed with cells containing 512 atoms for 10 ps, which ensured convergence of the free energy difference 
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arm to within 1 meV/ atom. In the second step we calculated i^vib~-^ref where, since the fluctuations 
in the energy differences U — C/j-ef were very small, we were able to use the second order formula (Eq. P). 

The problem in the calculation of thermal averages for a nearly harmonic system is that of ergodicity. 
For an harmonic system different degrees of freedom do not exchange energy, so in a system which is close 
to being harmonic the exploration of phase space using molecular dynamics can be a very slow process. 
We solved this problem following Ref. ||l^ whereby the statistical sampling was performed using Andersen 
molecular dynamics |^5[, in which the atomic velocities are periodically randomised by drawing them from 
a Maxwellian distribution. This type of simulation generates the canonical ensemble and overcomes the 
ergodicity problem. 

All the calculations were performed on a 64-atom cell with kpoints in a 7 x 7 x 7 grid for the high P/T 
state points and a 9 x 9 x 9 grid for the low P/T state points equivalent to 172 or 365 points in the IBZ 
respectively. 

The anharmonic contribution to the free energy of the solid turns out to be very small, being positive 
and equal to only a few meV/atom at low pressure and approximatively -20 meV/atom at high pressure. 



4.4 Error estimates for Fgoi 

The errors in Fperf are less than 1 meV/atom, the errors in -Fharm are ~ 3(4) meV/atom at low(high) 
P/T and the errors in Fahnarm are ^ 1(4) meV/atom at low(high) P/T; the total errors in Fgoi are « 3(6) 
meV/atom at low(high) P/T. 



5 Results and discussion 

We display in Figure 1 our calculated melting curve compared with the experimental zero pressure 
value 0], the DAC high pressure results fl^ , ^ and the high pressure shock datum [|l^. We also 
report in Figures 2a, 2b and 2c the volume change on melting, Vm, the entropy change on melting, Sm, 
and the melting gradient, dTm/dP, respectively. The errors in the melting curve arise from the errors in 
the calculated free energies and are ~ 50(100) K in the low(high) pressure part of the diagram respectively. 

The overall agreement with the experiments is extremely good; however, the low pressure results differ 
by more than 15 % (at zero pressure, 786 K compared with the experimental value of 933 K). Indeed, 
at zero pressure the agreement between the calculated and experimental volume change on melting and 
dTra/dP is rather poor (see Table ||). In addition, our calculations are not in very good agreement with 
the previous calculations of de Wijs et al. [^], although this is not necessarily surprising, since these latter 
calculations were based on LDA, while ours are based on GGA. Nevertheless, one might expect the results 
from LDA and GGA to be similar, since Al is a nearly free-electron like metal and therefore one would 
expect a very good DFT description with both LDA and GGA. To explore a possible reason why GGA 
does not predict the melting properties of aluminium very accurately we consider the zero pressure crystal 
equilibrium volume. This is predicted by GGA to be ~ 2% larger than the experimental value; this means 
that the calculated pressure for the experimental zero pressure volume is ~ -1-1.6 GPa. To see how this 
error propagates in melting properties we may devise a correction to the Helmholtz free energy such that 
the pressure is rectified: 

Fcorr = F + 6PV, (37) 

with 5P = 1.6 GPa. Using -Fcorr in our calculations we found the corrected melting curve, represented by 
the dotted line in Fig. 1, where we assumed 5P to be the same in the whole P/T range. The zero pressure 
corrected melting temperature is 912 K, which is in very good agreement with the experimental value 933 
K. The corrected volume change on melting, entropy change on melting and dT^/dP are also in much 
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better agreement with the experimental numbers. The correction is less important at high pressure, where 
dTm/dP is smaller. 

This point may be further illustrated by looking at the zero pressure phonon dispersion curves for 
Al. Since phonon frequencies depend on the interatomic forces, their correctness is surely important in 
the context of melting. In Figure 3 we display the GGA calculated phonon dispersion curves compared 
with experimental data |46|. Our calculations were performed both at the GGA zero pressure equilibrium 
volume and the experimental volume (both at 80 K). We notice that the agreement is good (though not 
perfect) if the calculations are performed at the experimental volume, and rather poor if the calculated 
zero pressure GGA volume is used instead. This indicates that GGA will probably yield better results if 
the GGA pressure is corrected in order to match the experimental data. 

In their work, de Wijs et al. [^] found good agreement between LDA and experiments. In their case 
a corrected LDA would lower the zero pressure melting point below 800 K. In order to understand this 
apparent different behaviour between LDA and GGA we have also calculated the phonons using LDA 
at the calculated equilibrium volume and also at the experimental volume (both at 80K). These are also 
reported in Fig. 3. In accord with previous LDA calculations [47| we found very good agreement with the 
experiments when the phonons are calculated at the LDA zero pressure volume, but the agreement becomes 
poor at the experimental volume, which is consistent with the result for the melting temperature |19]. 

In conclusion, both GGA and LDA predict an incorrect equilibrium volume at a fixed pressure, although 
LDA yields very good results for both the phonon dispersion curves and the zero pressure melting properties 
(which is probably accidental). For GGA the incorrect equilibrium volume propagates to an incorrect 
description of the phonon frequencies and the melting properties. If the GGA pressures are corrected so as 
to match the experimental data, the phonon dispersion and the melting properties come out in very good 
agreement with the experiments. These two behaviours are internally consistent, but point to an intrinsic 
error due to the use of GGA. Quantum Monte-Carlo (QMC) techniques [48| have been shown to predict 
the energetics with much higher accuracy than DFT [^], and calculations for systems containing more 
than 100 atoms have already been reported [^0|. We believe that in the near future it will be possible to 
use QMC for more accurate calculations of free energies. 

To summarise, we have calculated the melting curve of aluminium entirely from first-principles within 
the DFT-GGA framework. Our work is based on the calculation of the Gibbs free energy of liquid and 
solid Al, and for each fixed pressure the melting temperature is determined by the point at which the 
two free energies cross. Our results are in good agreement with the available experimental data, although 
they reveal an intrinsic DFT-GGA error which is responsible for an error of ~ 150K in the low pressure 
melting curve. This error is probably due to the incorrectly predicted pressure by GGA, and it becomes 
less important in the high pressure region, as dT^/dP becomes smaller. 
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Experiment 


LDA 


GGA 


GGA - corrected 


Tm (K) 


933 


890 (20) 


786 (50) 


912 (50) 


Sfn (^b) 


1.38 


1.36 (4) 


1.35 (6) 


1.37 (6) 


K^ (A3) 


1.24 


1.26 (20) 


1.51 (10) 


1.35 (10) 


dTm/dP (K GPa-i) 


65 


67 (12) 


81 


71 



Table 1: Comparison of ab initio and experimental melting properties of Al at zero pressm'e. Values are 
given for the melting temperature, Tm, entropy change on melting, Sm, volume change on melting, Vm, 
and melting gradient dTm/dP. The LDA results are from Ref. [19|; the experimental values for T^, Sm 



and dTm/dP are from Refs. |11], [52| and |51| respectively, and the experimental melting volume, Vm 
is calculated using the Clapeyron relation, Vm = SmdTm/dP. 
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Figure 1: Comparison of melting curve of Al from present calculations with previous experimental results. 
Solid curve: present work; dotted curve: present work with pressure correction (see text); diamonds and 
triangles: DAC measurements of Refs. [12] and [^] respectively; square: shock experiments of Ref . [p!^]. 
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Figure 2: Calculated pressure dependence of the melting properties of Al: a) volume change on melting, 
b) entropy change on melting and c) melting gradient. Solid curve: present work; dotted curve: present 
work with pressure correction (see text). 
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Figure 3: Comparison of the phonon dispersion curve for Al from present calculations with previous 
experimental results. Solid curves: present work with GGA; dotted curves: present work with GGA and 
with pressure correction (see text); dashed curves: present work with LDA; dot-dashed curves: present 
work with LDA and with pressure correction (see text); diamonds: experiments from Ref. [^]. 
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